Hydrodynamic Model for Particle Size Segregation in Granular 

Media 



Leonardo Trujillo 1 and Hans J. Herrmann 1 ' 2 
P.M.M.H.* , Ecole Superieure de Physique et de Chimie Industrielles, 10, rue Vauquelin, 75231 

Paris Cedex 05, France 
2 ICA-1, University of Stuttgart, Pfaffenwaldring 27 D-70569 Stuttgart, Germany 

(February 1, 2008) 

We present a hydrodynamic theoretical model for "Brazil nut" size segre- 
gation in granular materials. We give analytical solutions for the rise velocity 
of a large intruder particle immersed in a medium of monodisperse fluidized 
small particles. We propose a new mechanism for this particle size-segregation 
due to buoyant forces caused by density variations which come from differ- 
ences in the local "granular temperature". The mobility of the particles is 
modified by the energy dissipation due to inelastic collisions and this leads to 
a different behavior from what one would expect for an elastic system. Using 
our model we can explain the size ratio dependence of the upward velocity. 

PACS numbers: 45.70.Mg, 05.20.Dd, 05.70.Ln 



I. INTRODUCTION 

The physics of granular materials is a subject of current interest [|I]]. A granular medium 
is a system of many macroscopic heterogeneous particles with dissipative interactions. One 
of the outstanding problems is the so-called "Brazil Nut effect" 0: When a large intruder 
particle placed at the bottom of a vibrated bed tends to the top. This size segregation is due 
to the nonequilibrium, dissipative nature of granular media. Granular materials are handled 
in many industries. Many industrial machines that transport granular materials use vertical 
vibration to fluidize the material, and the quality of many products is affected by segregation. 
Size segregation is one of the most intriguing phenomena found in granular physics. A 
deeper understanding of this effect is therefore interesting for practical applications, and 
also represent theoretical challenge. 

A series of experiments |3|-|9| and computer simulations P,^0|-|17|] have elucidated differ- 
ent size segregation mechanisms, including vibration frequency and amplitude |^-^,TD,I3 



particle size P,^,|6|, p!0| , pd| , p!3| , p!4| , p!7| and size distribution [p!2|,|15| ; particle shape |16j; and other 
properties such as density ||^,|T3|,|T7]] and elastic modulus JT5 . 



Several possible mechanisms for size segregation have been proposed. One is segregation 
in presence of convection observed experimentally in three dimensions by Knight et ai, 
Q, and by Duran al, in two dimensions |J, under conditions of low amplitude and high 
acceleration vibration. In this case, both intruder and the small particles are driven up along 
the middle of the cell, and while the smaller particles are carried down in a convection roll 
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near the walls the intruder remains trapped on the top. In experiments performed by Vanel 
et al., in three dimensions, they observed two convective regimes separated by a critical 
frequency M. The first regime is associated with heaping and the second regime is similar 
to the one observed in Ref . [|J . Also, they reported a nonconvective regime observing a size 
dependent rise velocity. Employing large molecular dynamics simulations in two dimensions 
Poschel and Herrmann 0], and in three dimension Gallas et al, [15|, have recovered several 



aspects that are seen in experiments and recognize the lack of a theoretical description of 
the exact mechanism driving the segregation and the role of convection. 

Other segregation mechanism is associated to the percolation of small grains. Based on 
a Monte Carlo computer simulation, Rosato et al, |2| argue that each cycle of the applied 
vibration causes all the grains to detach from the base of the container. Then, the smaller 
particles fall relatively freely, while the larger particles require larger voids to fall downwards. 
The large grains therefore effectively rise through the bed. In the context of large-amplitude, 
low-frequency vertical shaking process (tapping), Jullien et al, predict a critical size ratio 
below which segregation does not occur pT0| , pLT|1 . This provoked some controversy [pj|-|2~I| and 
this threshold may be an artifact of the simulation model based on the "steepest descent 
algorithm" ||10|| . Experiments in Helle Shaw cells @J^,§] observed an intruder size dependent 
behavior, where the segregation rate increases with the size ratio between the intruder and 
the surrounding particles. Duran et al, formulated a geometrical theory for segregation 
based on the arching effect ||. They also claim experimental evidence for a segregation size 
threshold |J. In this picture the intruder contributes to the formation of an arch sustained 
on small grains on both sides. Between each agitation the small particles tend to fill the 
region below the arch. So, at each cycle the small particles move downward and the intruder 
effectively rises. Using a modification of the algorithm proposed by Rosato et al, Dippel 
and Luding find a good qualitative agreement with the non-convective and size-dependent 
rising 
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In another context, Caglioti et al, considered the geometrical properties of mixtures in 



the presence of compaction ||16|| . They established a relation for the effective mobilities of 
different particles in heterogeneous situations. 

The effect of the intruder density was studied by Shinbrot and Muzzio || . They observed 
an oscillating motion of the intruder on the top, which corresponds to the "whale effect" 
predicted by Poschel and Herrmann JTI3] . Also, they observed a reverse buoyancy in shaken 
granular beds. Ohtsuki et al, performed molecular dynamics simulations in two dimensions 
and studied the effects of intruder size and density on the height, and found no segregation 
threshold flTj|]. Recently Shishoda and Wassgren performed two dimensional simulations to 



model segregation in vibrofluidized beds []TT[ . They reported an height dependence with 
the density ratio between the intruder and the surrounding particles. In their model the 
intruder position results from a balance between the granular pressure (buoyant force) within 
the bed and the intruder weight. Their approach is in some sense similar to the model that 
we propose in this article. 

Subject to an external force, granular materials locally perform random motions as a 
result of collisions between grains, much like the molecules in a gas. This picture has in- 
spired several authors to use kinetic theories to derive continuum equations for the granular 
flow-field variables p2|-[29|j. Some of these theories have been generalized to multicompo- 
nent mixtures of grains PDHSHI . For different size particles in the presence of a temperature 
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gradient, Arnarson and Willits, found that larger, denser particles tend to be more concen- 



trated in cooler regions [0. This result was confirmed by numerical simulations |34|,|35 



However, this mechanism of segregation is a natural consequence of the imposed gradient of 



temperature and its not related to the nature of the grains |35[ . 

In this article we address the problem of size segregation using a kinetic theory approach 
in two and three dimensions (D = 2,3). We consider the case of an intruder particle 
immersed in a granular bed. We propose a segregation mechanism based on the difference 
of densities between different regions of the system, which give origin to a buoyant force 
that acts on the intruder. The difference of densities is caused by the difference between 
the mean kinetic energy among the region around the intruder and the medium without 
intruder. The dissipative nature of the collisions between the particles of a granular media 
is responsible for this mean energy difference, and modifies the mobility of the particles. 

The plan of this article is as follows. In Sec. II we derive a continuum formulation for 
the granular fluid, and introduce the definition of the "granular temperature" . In Sec. Ill 
we propose an analytic method to estimate the local temperature in the system. In Sec. IV 
we introduce the coefficient of thermal expansion. In Sec. V explicit solutions of the time 
dependence of height and velocity of the large particle are calculated. We can explain the 
size ratio dependence of the rise velocity and address the issue of the critical size ratio to 
segregation. To validate our arguments we make comparisons with previous experimental 
data. 



II. CONTINUUM FORMULATION 

We consider an intruder particle of mass mi and radius 77 immersed in a granular bed. 
The granular bed is formed of N monodisperse particles of mass mp and radius r>. The 
particles are modeled by inelastic hard disks [D = 2) or spheres (D = 3) in a /^-dimensional 
volume V = L D of size L. The size ratio is denoted = rijrp. The particles interact via 
binary encounters. The inelasticity is specified by a restitution coefficient e < 1. We assume 
this restitution coefficient to be a constant, independent on the impact velocity and the 
same for the fluid particles and the intruder. The post collisional velocities v' are given in 
terms of the pre-collisional velocities v by 

Vi,2 = Vi, 2 T [(vi - v 2 ) • n]n, (1) 

where the labels 1 or 2 specify the particle, n is the unit vector normal to the tangential 
contact plane pointing from 1 to 2 at the contact time, and the reduced mass m re d = 
m\m<il (mi + m 2 ). To calculate the dissipated energy we consider that energy is dissipated 
only by collisions between pairs of grains. In a binary collision the energy dissipated is 
proportional to AE = —m re d(l — e 2 )v 2 /2, where v is the mean velocity of the particles. 

In this work we use a generalized notion of temperature. In a vibrofluidized granular 
material a "granular temperature" T g can be defined to describe the random motion of the 
grains and is the responsible for the pressure, and the transport of momentum and energy 
in the system p6| . The granular temperature T g is defined proportional to the mean kinetic 



energy E associated to the velocity of each particle 
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- T s = ^ = ^LU m ^ • ( 2 ) 



2 y N N^\2 

We expect a continuum limit to hold for N ^> 1, when the small particles may be 
considered as forming a granular fluid. In order to develop an analytic study, we assume 
that the uniformly heated granular fluid can be described by the standard hydrodynamic 
equations derived from kinetic theories for granular systems ||22|| . In this study, we focus on 
a steady state with no macroscopic flow. 

The balance equation for the energy is 

V • q = - 7 , (3) 

where q is the flux of energy and 7 is the average rate of dissipated energy due to the 
inelastic nature of the particles collisions. The constitutive relation for the flux of energy, 

q = -KVT g , (4) 

defines the thermal conductivity k. Consequently, we have 

V • (kVT 9 ) = 7. (5) 

A uniformly fluidized state can be realized when the granular material is vibrated in the 
vertical direction, typically as z (t) = A sin(uj t), with the amplitude A and the frequency 
ujq = 2nf, so that one can define a typical velocity u = A q ujq. In the experiments the 
excitation is described by the dimensionless amplitude T = A u /g, where g is the gravi- 
tational acceleration. As a first approximation the effect of the external force experienced 
by the particles due to the gravitational field is neglected in the description of the granular 
flow. Experimentally this corresponds to the regime T ^> 1. So, the momentum balance, 
in the steady state, implies that the pressure p is constant throughout the system. 

The hydrodynamic equations close with the state equation, the collisional dissipation 7 
and the transport coefficients for a granular medium. In the limit N ^> 1 the constitutive 
relations are determined as function of the properties of the small grains. The transport 
coefficients are assumed to be given by the Enskog theory for dense gases in the limit of 
small inelasticity. 

The total pressure should be essentially equal to that of the small particles, the contri- 
bution of the intruder being negligible, since N 3> 1- For a dense system the pressure is 
related to the density by the virial equation of state, which in the case of inelastic particles 
is 



P = ~^ nT 9 



1 + ^n 9o (2r F ) D 



(6) 



where n = N/V is the number density of small grains, tt D = 2n D / 2 /T(D/2) is the surface 
area of a D-dimensional unit sphere, v is the area (volume) fraction v = Qj^nrp /D, and go 
is the pair correlation function for two fluid particles. In two dimensions the pair correlation 
function is [55 1 



- v) 
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with the area fraction v = nnrp. In three dimensions the pair correlation function is [p7 



(2-f) 

with the volume fraction z/ = 47rnr F /3. 

The state-dependent thermal conductivity possesses the general form ]22| 



k = KoyTg, (9) 

where the prefactor k is a function of the fluid particle properties, and can be calculated 
using a Chapman-Enskog procedure through the solution of Enskog transport equation 
23| , |2T|p9| , |3"8"| . The explicit expressions of these prefactors are given in Appendix A. 



To estimate the collisional dissipation rate 7 we consider the loss of average kinetic energy 
per collision and per unit time. In a binary collision the kinetic energy dissipated can be 
expressed in terms of the granular temperature as AT g = —(1 — e 2 )T g /2. For the fluid 
particles, the average collision frequency up is proportional to iOp ~ an d we assume 

that it is given by the Enskog collision frequency [[$| 

^-^^-(AyV. do) 



This form for the frequency of collisions is justified for a granular medium. This is a conse- 
quence that the average spacing between nearest neighbor s is supposed to be less than the 
grain diameter (s <C 2r>) ||22|| . Multiplying AT g by the collision rate up and the number 



density n = N/V, we obtain the collisional dissipation rate 7^ for the fluid particles 

O / 9 \ V 2 



^=(1 - e 2 )n 2 g (2r F ) D - 1 (-=-) T^. (11) 



If = , — , , . , 

2V27T \m F / 

In order to simplify the mathematical notation let us express jp as 

If = ^Tf 2 , (12) 

where the dissipation factor £p contains the prefactors which multiply T^ 2 in Eq . (|IT]) , this 
is 

^-| ( i- e >^(2. f rW 1/2 (13) 

Zy Zn \mp / 

To understand the essential features of the intruder's presence in the granular medium, 
it is adequate to adopt a simplified point of view. If the mean velocity of the fluid particles 
is u, the flux of fluid particles which strikes the intruder's surface can be estimate as nu. 
Multiplying this flux by the area of the intruder Vtp,^ , we can calculate the number of 
fluid particles which strike the surface of the intruder per unit time, and written in terms of 
the granular temperature we have 

tt D d-1 f m i + m F\ 1/2 rpl/2 /-mN 
V27T V mimp J y 



5 



So, the local density of kinetic energy dissipated in the region near the intruder is 

'i-e^ g „( rF + n)-'f^±^V /2 rr. US) 



77 2a/2tT^ -/yauvr -x, y mi m F J ~ 9 

In the simplified form Eq.(|l5|) can be expressed as 

li = tiT^ 2 , (16) 
where the dissipation factor £/ is defined as 

«'=dfe (1 - e V°^ + ^ fer) ■ (17) 



III. LOCAL TEMPERATURE DIFFERENCE 

The intruder's presence modifies the local temperature of the system due to the collisions 
that happen at its surface. The number of collisions on the surface increases with the size 
of the particle, but the local density of dissipated energy diminishes. From Eq.(|5]) we can 
calculate within a sphere of radius tq the value of the temperature in the granular fluid in 
presence of the intruder and compare it with the temperature in the granular fluid without 
intruder, we will denote these temperatures T\ and T 2 respectively (see Fig.(l)). This is a 
simple method to estimate the temperature difference between a region with intruder and a 
region without intruder AT g = T\ — T 2 . 



. OO^ 
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(a) 



(b) 



FIG. 1. Schematic representation of the regions used to calculate the granular temperature, (a) 
Region around the intruder within a sphere of radius r$ and (b) region without intruder. 

Let us concentrate on solutions with radial symmetry. The solutions of Eq.(^|), for an 
arbitrary dimension D, satisfy the equation 



1 



dr 



£ T 3/2_ 



(18) 



This nonlinear differential equation can be simplified by the fact that the pressure is consid- 



ered constant throughout the system and remembering that p 
the resulting equation may be written in terms of w = T^ 2 , 



T„. So, linearizing Eq. 
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dr 2 r dr 

where 



d 2 w D — 1 dw , 9 

+ — — = A V (19) 



The collisional dissipation rate can be decomposed in two parts. We propose this de- 
composition supposing that the energy dissipation around the intruder is dominated by the 
collisions between the small grains and the intruder, then the dissipation rate in this region 
is given by Eq.flTJ). In the rest of the system the dissipation rate is dominated by the 
collisions between small grains only. In this case the dissipation is given by Eq.flTTl). 

First, let us consider the "inhomogeneous case" when the intruder is localized in the 
center of the system (r = 0), see Fig. (la). The dissipation factor £ can be decomposed 
in two parts: £ = £/ for the region near the intruder (r = 77), and £ = £_p for the region 
(rj<r< To), where r$ is the radius of the considered region. 

For the inhomogeneous case we express Eq.flTPP as 

,0 n 1 j f X 2 T w for < r < rr, 

d z w D - 1 dw \ 1 - ' , . 

*? + — L . , < < 21) 

^ \ F w for rj < r < ro, 

where, Xj = £//2k and = £i?/2/«o- The solution of Eq. (pl|) is determined by the boundary 
conditions imposed upon the system. As boundary condition we suppose that the system 
is enclosed by an external surface of radius r at temperature T g (r ) = T Q (respectively, 
w(r ) = w ). 

Let us denote Ti_(r) the granular temperature for the region (0 < r < 77), and Ti + (r) 
the granular temperature for the region (77 < r < ro) (respectively, W-(r) = JTi_{r) and 



u>+(r) = wT 1+ (r)). The intruder's presence imposes internal boundary conditions. On the 
inner surface, the temperature should satisfy 

™-(r)\ r=ri = w + (r)\ r=ri . (22) 

The flux of energy also imposes another internal boundary condition. If we suppose the 
flux of energy continuous on the inner surface, from Eq.(^) the granular temperature should 
satisfy 

dw_(r) _ dw+(r) 

dr dr 

r=rj r=rj 



A. Solution for 2D 



The solutions to Eq.(|2lD for D = 2 are a linear combination of the modified Bessel 
function of order zero W\(r) = {Io(Xr), K (\r)}. The general solution is 

w_(r) = A_J (A 7 r) + B^K (Xjr) for < r < n, (24) 
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and 



w + (r) = A + I (X F r) + B + K (X F r) for 77 < r < r , 

where v4_, A+, £?_ and 5 + are constants that must be determined from the 
conditions. 

The function K (Xr) diverges when r — * 0, then 

= 0. 

When r = r the Eq.(p5|) should satisfy the boundary condition 



this is, 



w + (r)\ = w , 



A + I (X F r ) + B + K (X F r ) = w . 



(25) 
boundary 

(26) 
(27) 
(28) 



On the inner surface the boundary condition fl22|) W-(rj) = w + (rj) leads to 



A_/ (A / r / ) = A + I Q (X F n) + B + K (X F rj) 



/o(Ajrj) 

The inner boundary condition (|23"D leads to 



/o(A/rj) ' 



A-Xih(Xin) = A + X F Ii(X F rj) - B+XpK^Xprj) 



A. 



Equating Eqs.(g9|) and (|3Cf) we find 



A h(X F rj) B K^Xprj) 



hi.hr 1) 



A+ = A F / (A J r J )ir 1 (A J ,r J ) + A / / 1 (A J r J )K (A F r J ) 
B + AF/o(A/r/)/i(A F r/) - A//i(A7-r/)/ (AFr/) 
= Qab 



From Eqs. (|28]) and ([H]) the constant B + should be 

w 



B 4 



&ABlo(^Fro) + K (X F r ) 
Substituting Eq.(|32|) into fl31~|) we have 

A _ WqQab 

+ <S>ABlo(>^Fro) + K {X F r Q ) 

Substituting Eqs.(B2) and (B3) into (E3) we have 



A. 



w 



QABlo(^Fro) + K (X F r ) 



6 



AB 



I (Xiri) /o(A/r 7 ) 



(29) 



(30) 



(31) 



(32) 



(33) 



(34) 



The granular temperature in the inhomogeneous case is 



Ti(r) 



(A_/ (A/r)) 2 for < r < 77, 

(A + I (X F r) + B + K (X F r)) 2 for 77 < r < r , 



(35) 



where the constant A_, A + and 5+ are given by the Eqs.(|J), (q3|) and ( J32|) , respectively. 

In the "homogeneous case", see Fig. (lb), the prefactor Xj = 0. Then the granular 
temperature T 2 (r) is 



T,(r) 



' Io{XFr) ) T v . 



(36) 



Now we are interested in determining the temperature difference AT g between case 1 
and 2 in the granular fluid. For this we calculate the granular temperatures at r = 0. When 
r — > the modified Bessel function of zero order tends to 1 f39 |. So, in Eq.(^) and (|36|) in 
the limit r->0we find that 



(Jo^r)) 2 ~ 1, 
(I (X F r)) 2 ~ 1. 

Then, the temperature difference is 

1 



(37) 



QABlo(^Fro) + K (X F r ) 



IoQfTj) + K o(X F n) 

Jo(A z rj) /o(A7-r 7 ) 



To 

(38) 



in two dimensions. 



B. Solution for 3D 

When D = 3, the solution of Eq.(|2lD is given in terms of the spherical modified Bessel 
functions of zero order wi(r) = {io(Xr) = sinh(Ar)/Ar, ko(Xr) = e~ Xr /Xr}. The general 
solution in this case is 



u>_(r) = A_i (A/r) + S_fc (A/r) for < r < 77, 



and 



w+(r) = A + i (X F r) + B + k (X F r) for 77 < r < r . 
The function fco(Ar) diverges when r — > 0, then 

= 0. 



(39) 



(40) 



(41) 



The constants A_, A + and 5 + , are calculated from the boundary conditions in a similar 
way as before. 
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A. 
A, 



w 



®ABio(^Fr Q ) + k (X F r ) 

WQ&AB 

@ab«o(Af»"o) + k (\ F r ) 
w 



i (X F ri) k (X F ri) 
VAB—p: r + —77 r 



©ab^A^o) + /^(Ai^o)' 
where in this case the factor Qab is 

A F z (A / r / )/ci(A F r / ) + A / zi(A / r / )/c (AF^/y 



(42) 

(43) 
(44) 



e 



AB 



AF«o(A/r/)ii(A F r/) - A/i 1 (Ajr/)io(A i ?r/) ' 
The granular temperature in the inhomogeneous case in 3D is 

(A_i (A 7 r)) 2 for < r < 77, 

(A + i (A F r) + B + k (X F r)) 2 for r 7 < r < r , 



(45) 



TAr) 



(46) 



where the constant A_, A + and £>+ are given by the Eqs.([42|), ([43|) and 

In the "homogeneous case" the prefactor A/ = 0. Then the granular temperature T 2 (r) 

is 



T,(r) 



v «o(Airr ) / 



(47) 



Again the temperature difference AT 9 is calculate at r = between case 1 and 2. When 
r — > the spherical Bessel function of zero order tends to 1 |39|, then 



/ sinh(A/ 



r) 



sinh(Aj? r) 



(48) 



Then, the temperature difference is 

a:/;, ' 1 



,©ab«o(Af^o) + k (X F r ) 
in three dimensions. 



6 



io(^Fri) , k {X F ri) 



1 \ 2 



?o(A/r/) i (A/r 7 ) 



^o(Af^o) 



To (49) 



C. Energy equipartition breakdown 

Let us define the temperature ratio r = Ti(0)/T 2 (0). In two dimension we have 



r 



/o(A F r ) 



s 6ab/o(Af7"o) + K (Ajrr ) 
and for three dimensions, 



Q Io(A F rj) KojXpn) 
AB J (A/r/) / Q (Ajr/) 



1 \ 2 



(50) 
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_( ^o(Af^o) 

y&ABioi^Fro) + ko(X F r ) 

since X F > A/ we can verify that Ti(0) > T 2 (0), this means r > 1. So, the temperatures 
ratio between the region with intruder and the region without intruder are different. In our 
model this lack of equipartition is due to a difference between the collisional dissipation rate 
related to the particle sizes. In the elastic limit e — ► 1 the energy equipartition is restored 
r — ► 1. In Fig. 2, we present the qualitative behavior of r with the size ratio <p — r i/ r F, 
for different values of the coefficient e. The granular temperature difference increases with 
and depends on e. We can see that r is nearly constant and very close to unity when 
e = 0.99. 



i (X F ri) k (X F rj) 

~* t\ — \" + ~~7\ — r 

^o(A J rr/) zo(Airj) 



(51) 




Recently, this quantity was directly measured in experiments performed by Feitosa and 
Menon POfl . They observed that energy equipartition does not generally hold for a binary 
vibrated granular system. They reported that the ratio of granular temperatures depends on 
the ratio of particle mass densities. Also in fluidized binary granular mixtures the breakdown 
of the energy equipartition was observed experimentally |4l| and described theoretically in 
the framework of the kinetic theory E2| 



IV. THERMAL EXPANSION 



The definition of the granular temperature allows us to establish a statistical mechanics 
description of the granular medium. The change of mean energy of the system is basically 
due to a mechanical interaction with their external parameters (e.g., the amplitude A Q and 
the frequency Uq = 2irf of vibration, the volume of the system V, and the pressure p). 
Theoretically we can derive the energy relaxation to the steady state for a driven granular 
medium ||43|| . The system increases its energy as a result of external driving while its 
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decreases its energy by dissipation. The work W done to change the volume of the system 
from V to a certain quantity V + dV is equal to the change of its mean energy and its related 
to the mean pressure and volume by dW = pdV + Vdp. From the definition of granular 
temperature, the change of the granular temperature depends on the mean kinetic energy of 
the particles. A volume change dV is related to a temperature change dT g by the equation 
of state (|). 

We can express V as a function of T g and p, V — V(T g ,p). Thus given infinitesimal 
changes in T g and p, we can write 

dV =(m)/ T ° + (w)r, dp ' (52) 
= aVdT g — kpVdp, 

where a is the thermal expansion coefficient defined as 

1 ( dV\ 1 ( dn\ 



and k p is the "isothermal compressibility" defined as 

(54) 

If in a first approximation we neglect the variations of the coefficients a and k p , we can 
integrate Eq.fl5"2"D and find 




V{T g , p) = V exp [aAT g - k p Ap] , 
^V [l + aAT g - k p Ap] . 

From the temperature difference, Eqs.(|38D and fl49|) , we can conclude that we also have a 
change in the effective density (p = nmj) of the granular fluid between case 1 and 2, through 
the thermal expansion of volume produced by the difference of temperatures, Eq. (|55D , thus 



p = p(l-aAT g ). (56) 

The thermal expansion coefficient can be derived from the equation of state (|]) and 
definition The general form of the coefficient of thermal expansion is 

a = ic?(i/), (57) 

where C(v) is a correction due to the density of the system. In the dilute limit v —>■ and 
C{y) — > 1, and the above expression tends to the expected value for a classical gas a = l/T Q . 
The explicit form of C{v) is given in Appendix B. 
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V. SEGREGATION FORCES 



Now we propose that this density difference leads to an effective buoyancy force f^, similar 
to the Archimedean force 

f 6 = ApV lg , (58) 

where Ap = —apAT g , Vj = ^§-rf is the D-dimensional volume of the intruder and g is the 
gravity field. The intruder also experiences a viscous drag of the granular fluid. The drag 
force id is considered to be linear in the velocity of segregation u(t), and is like the Stokes' 
drag force 

f d = -&Kfir T u(t), (59) 

where /i is the coefficient of viscosity of the granular fluid The state-dependent viscosity 
posesses the general form [p2| 



H = HoyTg, (60) 

where the prefactor /i is a function of the fluid particle properties, and can be calculated 
using a Chapman-Enskog procedure for the solution of Enskog transport equation. The 
explicit expressions of these prefactors are given in Appendix A. 

Equations (|58|) and ( |59|) express the acting forces in the segregation process 



f seg = f b + U (61) 
Therefore, the equation of motion that governs the segregation process is 
Q D D du{t) VL D D 

~D ri P ~dT = ~~D ri apAT 9 g ~ 67r ^/u(t). (62) 

Now we suppose the granular system contained between two large parallel plates per- 
pendicular to the gravitational field. We take the reference frame positive in the upward 
vertical direction. Arranging terms in Eq.(^) we find the following differential equation 

du{t) QTrDp{(t)r F ) l - D 

-^r = aAT ° 9 n^P — (63) 

where we have expressed the intruder's radius as function of the size ratio dependence 
rj = 4>rp, and the solution of this differential equation is the rise velocity of the intruder 



aAT g gt 

where the time-scale to is given by 



exp 



a-D * 



(64) 



to s eJ^p- (65) 

The force balance between the drag force and the buoyant force gives the settling 
velocity u s 
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aAT q t 



a-D 



(66) 



The time dependent intruder height z{t) is 

aAT g gt 



z(t) 



t-t n 



exp 



■l-D 



tr 



(67) 



On a qualitative level our model satisfactorily reproduces the observed phenomenology: 
a large intruder migrates to the top of a vibrated bed, and the rise velocity increases with the 
intruder size. The solutions ( |64"D and flBTD are plotted in Figs. 3 and 4. Our results resembles 
the experimental intruder height time evolution described in Refs. @ and 0. However 
the model can not describe the intermittent ascent of the intruder since we calculate the 
mean velocities. Using the following model parameters: mass particle density of 2.7 gem -3 
(Aluminum), r F = 0.1 cm, e = 0.83, u = 0.34, N = 5 x 10 3 , g = 100 cms -2 , r = 2rj and 
T = 30 x 10 6 gcm 2 s~ 2 we obtain that the order of magnitude of z(t) (Fig. 4) coincides with 
the values reported by Cooke et al. (See Fig. 3 Ref. @). 
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FIG. 3. Intruder segregation velocity u(t). The parameters are: mass particle density of 2.7 
gem -3 (Aluminum), rp = 0.1 cm, e = 0.83, v = 0.34, iV = 5 x 10 3 , g = 100 cms -2 , tq = 2r/ and 
T = 30 x 10 6 gcm 2 s -2 . 
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FIG. 4. Intruder height time dependence z(t). The parameters are: mass particle density of 2.7 
gem -3 (Aluminum), rp = 0.1 cm, e = 0.83, v = 0.34, jV = 5 x 10 3 , g = 100 cms -2 , r = 2r/ and 
T = 30 x 10 6 gcm 2 s~ 2 . Inset: Measured intruder height (Fig. 3, Ref.[6]). 
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From the settling velocity u s ( |66"D we show explicitly the dependence on size. This is 
proportional to the size ratio <fi and the granular temperature differences AT g which also 
depends on the size ratio. Its agrees with the experimental fact that the larger the radius 
of the intruder, the faster is the ascent, reported by Duran et al. ||. The plotted solution 
(|66|), Fig. 5, fits well with the experimental values from Ref. || for <fi > 4, shown in Fig. 
5. In this experiment, Duran et al, claim the experimental evidence of a segregation size 
threshold at <p c = 3.3, below which the intruder does not exhibit any upwards motion. Our 
model's continuous aspect doesn't allow for the existence of this threshold. We argue that 
this discrepancy comes from the fact that experimental measures in this regime should be 
very difficult to earring out. They reported to not have observed any upward motion after 
one hour. If we see the height profile shown in Fig. 4 for = 2, we note that it is very close 
to zero after one minute and the slope is very low. Therefore this threshold should not exist. 

0.025 
0.02 

-55-0.015 

E 

-5" 0.01 
0.005 



" 2 4 6 8 , 10 12 14 16 

<N r , /r F 

FIG. 5. Intruder segregation velocity dependence on <j>. The parameters are: mass particle 
density of 2.7 gem -3 (Aluminum), rp = 0.75 cm, e = 0.85, v = 0.25, N = 5 x 10 3 , g = 100 cms -2 , 
ro = 2rj and To = 126 x 10 6 gcm 2 s~ 2 . The data points come from Ref. [5]. 

VI. CONCLUSIONS 

We derived a phenomenological continuum description for particle size segregation in 
granular media. We propose a buoyancy-driven segregation mechanism caused by the dis- 
sipative nature of the collisions between grains. The collisional dissipation rate naturally 
leads to a local temperature difference among the region around the intruder and the medium 
without intruder. In this model we proposed that the intruder's presence develops a tem- 
perature gradient in the system which gives origin to a difference of densities. The granular 
temperature difference is due to the fact that the number of collisions on the surface in- 
creases with the size of the intruder, but the local density of dissipated energy diminishes. 
So, the region around the intruder its hotter than the region without intruder. From this 
temperature difference we can conclude that we have a change in the effective density of 
the granular fluid. This leads to an effective buoyancy force that is the responsible for the 




16 



intruder's upward movement. 

In this work we made use of the tools of kinetic theory of gases to calculate the granular 
temperature. We observed a breakdown of the energy equipartition. And this is in agreement 
with other reported experiments and models. In a certain sense our theory unifies the 
different aspects observed in the size segregation phenomenon. Explicit solutions of the 
dependence of height and velocity are calculated. The geometrical effect of a segregation 
threshold is not supported by our model. The intruder size dependence appears naturally 
in our model. 

Very recently it has been shown experimentally [I3| and by computer simulations [J45[| , 
that the convection phenomenon in granular fluids comes from the effect of spontaneous 
granular temperature gradients, due to the dissipative nature of the collisions. This temper- 
ature gradient leads to a density variations. The convection rolls are caused by buoyancy 
effects initiated by enhanced dissipation at the walls and the tendency of the grains at the 
center to rise. So, this segregation mechanism could be described in the hydrodynamic 
framework proposed in this work subject to the appropriate boundary conditions. 
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APPENDIX A: TRANSPORT COEFFICIENTS 



In this appendix the prefactors appearing in Eqs.([|) and (|60D are derived. Using a 
Chapman-Enskog procedure for the solution of the Enskog transport equation, the transport 



coefficients for nearly elastic particles have been derived in Refs. p3| and |f2~7 |. 
In 2D the thermal conductivity k is [E3] 



k = 3nr F ( ] 

\m F J 



V2 r 11 3/ 16\_ 
1+ 3G + 4l 1 + d G 



T l/2 

a ' 



(Al) 



where G is ug , go is the 2D pair correlation function given in Eq.(^), and v is the area 
fraction v = nnr F . It is convenient to express Eq.( |Al| ) introducing the prefactor kq deffmed 

as 



k = 3nr F 



1/2 r 



m F J 



, 1 1 3 

1+ 3G + 4 



97T 



(A2) 



The result ( |A1[ ) takes the form 



K 



(A3) 



In 3D the thermal conductivity is p? 



15 / 7T \ 



1/2 



-nr F 



\m F J 



5 1 6 / 32\ „ 
2iG + A l + ^) G . 



T l/2 
± 9 ' 



(A4) 
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where G is is go, go is the 3D pair correlation function given in Eq.(H), and v is in this case 
the volume fraction v = Aimrp/3. In 3D the prefactor kq is defined as 



15 / 7T y/2 



5 1 6 / 32 \ ^ 



The shear viscosity /j in 2D is [E3 



(A5) 



/i = -nr^(vrmp) 1 / 2 



2 + - +(1 + -| G 



T 1/2 . 



It is convenient to express Eq. (IA3) introducing the prefactor /i defined as 



Aio = -nr F (7im F ) l/2 
So, the result (|A6|) takes the form 



2 + 



G 



7T 



G 



A 1 — ^o\Tg- 



In 3D the shear viscosity is [27] 



/j, = -nr F (7rm F ) 1 / 2 



T l/2 



and the prefactor /xq in 3D is defined as 

1 



yUo = -nr F (vrm F ) 1/2 



5 1 4 / 12\ „ 
1 + T6G + 5( 1 + V) G 



(A6) 



(A7) 



(A8) 



(A9) 



(A10) 



APPENDIX B: THERMAL EXPANSION COEFFICIENT 



We can consider the volume of the system as a function of the granular temperature and 
the pressure V = V(T g ,p). A change in the granular temperature dT g and the pressure dp, 
leads to the corresponding change in the volume dV 



<n = [^-) dT g +{^-\ dp. 



dp 



(Bl) 



As we have supposed that the pressuer of the system is more or less constant, we can 
approximate dp ~ 0. The increment of volume dV with an increment of the granular 
temperature dT g is 



(B2) 



Thus, 



dV ( dV 



dT g \dT g J p 



(B3) 
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or 



'dV_ 
dT„ 



dV 



(B4) 



and in terms of the number density n, we have 



dn 



dT, 



9/ p,N 



dn 



p,N 



(B5) 



From the definition of the coefficient of thermal expansion Eg. ([53"D , and from the above 



statement, we find 



1 f dn 

a = 



n \dT, 



9 J p,N 



dn 



p,N 



-1 



(B6) 



The partial derivative (dT g /dn) Pt N can be calculated from the equation of state (|B7|). In 
2D the equation of state is 



P = ^—^-nT g 



l + 2u 



[1-uf 



(B7) 



where v = nnr 2 F . So, an elementary calculation leads to 



'ST, 
dn 



2 p 8(z/ 3 -3z/ 2 -8z/-8)(z/- 1) 



From Eg. ( P36|) one obtains: 

2 n 

a 



is 2 + 



(B8) 



{v 2 + 8) 2 



(B9) 



1 + ep 8(z/ 3 - 3z/ 2 - 8u - 8)(u - 1) ' 
Using the equation of state ( |B7[ ) we can express a in function of the granular temperature 

1 (v 2 + 8) 2 



a 



T 9 (z/ 3 -3z/ 2 -8z/-8)(z/-l) 



this is 



1 



where the correction coefficient due to the density of the system is defined as 

CM = {u2 + 8)2 

y ' ~ (z/ 3 -3u 2 -8u-8){u- 1)' 
For three dimensions the equation of state is 



(BIO) 



(Bll) 



(B12) 
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1 + e 
2 



nT s 1 + 4i/ 



(2-^) 
2(1 - z/ 3 



(B13) 



In a similar way we find for 3D that the coefficient of thermal expansion is 



a = 



1 (z/ 3 -z/ 2 -z/- 1) 
7^(z/ 4 -4z/ 3 + 4z/ 2 + 4z/+l) 



(B14) 



and the correction coefficient C(v) in 3D is defined as 



C(u) 



(z/ 3 -z/ 2 -z/-l)(z/-l) 
(z/ 4 - 4z/ 3 + 4z/ 2 + 4z/ + 1) ' 



(B15) 



where z/ = 4n7rr 3 7 /3. 
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